// ///////////////////////////// MIT License //////////////////////////////////// //
//                                                                                //
// Copyright (c) 2010 David Zsolt Manrique                                        //
//                    david.zsolt.manrique@gmail.com                              //
//                                                                                //
// Permission is hereby granted, free of charge, to any person obtaining a copy   //
// of this software and associated documentation files (the "Software"), to deal  //
// in the Software without restriction, including without limitation the rights   //
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      //
// copies of the Software, and to permit persons to whom the Software is          //
// furnished to do so, subject to the following conditions:                       //
//                                                                                //
// The above copyright notice and this permission notice shall be included in     //
// all copies or substantial portions of the Software.                            //
//                                                                                //
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     //
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       //
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    //
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         //
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  //
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN      //
// THE SOFTWARE.                                                                  //
//                                                                                //
// ////////////////////////////////////////////////////////////////////////////// //

#include <iostream>
#include <ostream>
#include <sstream>
#include <vector>
#include <algorithm>
#include <complex>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <fstream>
#include <string>
#include <algorithm>
#include <iomanip>
#include <vector>
#include <list>
#include <set>
#include <cstdarg>
#include <stdexcept>

#include <Eigen/Eigen>
#include <cmdioutil.h>

#include "hamiltonian.h"

template<typename scattering>
void out_coef(double e, scattering & s,std::ostream & os,bool comment = false,bool advanced = false,int w = 18,int p = 10)
{
        //os.flags(std::ios::scientific);
        os.flags(std::ios::fixed);
        if(comment)
        {       
                std::ostringstream oss;  oss << "E-Ef(=" << std::setprecision(2)<< std::fixed << e <<")";
                os << '#';
                os << std::right << std::setw(w-1) << std::setprecision(p) << oss.str() << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Tran a->b" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Tran b->a" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Ref a->a" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "Ref b->b" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "a oc" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "b oc" << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << "max abs S*S'-I" << ' '; 
//                os << std::right << std::setw(w) << std::setprecision(p) << "eiva tt' a->b" << ' '; 
                os << std::endl;
        }
        else
        {   
    	    os.flags(std::ios::scientific);
            os << std::right << std::setw(w) << std::setprecision(p) << e << ' ';
            if(s.ano*s.bno !=0)
            {
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefba << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefab << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefaa << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << s.coefbb << ' ';
            } 
            else
            {
                os << std::right << std::setw(w) << std::setprecision(p) << 0.0 << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << 0.0 << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << (double)s.ano << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << (double)s.bno << ' ';
                s.checknorm = 0.0;
            }
            os << std::right << std::setw(w) << std::setprecision(p) << s.ano << ' '; 
            os << std::right << std::setw(w) << std::setprecision(p) << s.bno << ' '; 

            if(s.ano*s.bno ==0) s.checknorm = 0.0;
            os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << s.checknorm << ' '; 
            os << std::endl;                

	    os.flags(std::ios::fixed);

/*          for(int i = s.ttabe.rows()-1; i >= 0; i--)
            os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << s.ttabe(i).real() << ' ';  
            os << std::endl;                
            
            if(advanced)
            {
                os << std::endl;
                os << "advanced: scatterer wv for scatterer channles (experiment 2010jul16)" << std::endl;
                for(int i = s.ttabe.rows()-1; i >= 0; i--)
                os << std::right << std::setw(33) << std::setprecision(10) << std::scientific << s.ttabe(i).real() << ' ';  
                os << std::endl;                

                for(int i = 0; i < s.yml.rows(); i++)
                {
                    for(int j = s.yml.cols()-1; j >= 0; j--)
                    {
                        os << std::right << std::setw(15) << std::setprecision(5) << std::scientific << s.yml(i,j).real() << ' ';  
                        os << std::right << std::setw(17) << std::setprecision(5) << std::scientific << s.yml(i,j).imag() << ' ';  
                    }
                    os << std::endl;
                }
                os << std::endl;
                os << std::endl;
            }
*/

        }
        

}



template<typename lead,typename scatterer>
void get_amb(int argc, char* argv[],lead & a,scatterer & m,lead & b,double & emin_p, double & emax_p,double & ef_p,int & nume_p,double &tol_p,bool & sqrte_p,bool &advanced)
{
	std::cout << "# ";
	for(int i = 0; i < argc; i++)
	std::cout << argv[i] << " ";
	std::cout << std::endl;
	
        cmd::arguments arg(argc,argv);
        cmd::switcharg<false> help(arg,"-h","--help", "print help message");
        cmd::switcharg<false> info(arg,"-i","--info", "print info/example message");

        cmd::vararg<double> emin(arg,"-emin","--emin","min of e-range (ef will be set to zero), <arg>=-1.0",-1.0);
        cmd::vararg<double> emax(arg,"-emax","--emax","max of e-range (ef will be set to zero), <arg>=1.0",1.0);
        cmd::vararg<int> nume(arg,"-n","--nume","number of energy intervals, <arg>=1",1);
        cmd::vararg<double> efermi(arg,"-ef","--ef","set fermi energy explicitely",0.0);
        cmd::switcharg<false> sqrtesampling(arg,"-sqrt","--sqrt-sampling", "sqrt e point sampling around ef");
        
        cmd::vararg<util::pair_expression<int,int,':'> > ha(arg,"-ha","--ha","lead-a uc pair atom index range, 1:3 or -1:-4","");
        cmd::vararg<util::pair_expression<int,int,':'> > hb(arg,"-hb","--hb","lead-b uc pair atom index range, 1:4 or -1:30","");

        cmd::vararg<util::pair_expression<int,int,':'> > ah(arg,"-ah","--ah","lead-a uc pair atom index range (h taken from opposite side)","");
        cmd::vararg<util::pair_expression<int,int,':'> > bh(arg,"-bh","--bh","lead-b uc pair atom index range (h taken from opposite side)","");

//        cmd::varargxx<util::pair_expression<int,int,':'>,double > fa(arg,"-fa","--fa","fake lead-a layer atom index range and coupling, 1:2 10.0","",10.0);
//        cmd::varargxx<util::pair_expression<int,int,':'>,double > fb(arg,"-fb","--fb","fake lead-b layer atom index range and coupling, 1:2 10.0","",10.0);

        cmd::varargxx<util::pair_expression<int,int,':'>,std::string> la(arg,"-la","--lead-a","replace i:j h range and lead-a ham. file","","");
        cmd::varargxx<util::pair_expression<int,int,':'>,std::string> lb(arg,"-lb","--lead-b","replace i:j h range lead-b ham. file","","");

        cmd::vararg<double> tol(arg,"-tol","--tol","channel separation tolerance def.: 1e-8",1.0e-8);

        cmd::vararg<double> randa(arg,"-ra","--randa","randomize lead-a, to avoid singular coupling, <arg>=0.0 ",0.0);
        cmd::vararg<double> randb(arg,"-rb","--randb","randomize lead-b, to avoid singular coupling, <arg>=0.0",0.0);
        cmd::vararg<double> randab(arg,"-r","--randab","randomize both leads, to avoid singular coupling, <arg>=0.0",0.0);

        cmd::switcharg<false> homogenize(arg,"-ab","--ab", "homogenize the two leads(taking average)");        
        cmd::vararg<double> smallimag(arg,"-epsim","--epsim","add small imaginary in the scatterer diagonal ham. <arg>=+0.01",0.01);

        cmd::vararg<std::string> chain(arg,"-chain","--chain", "1D chain","");        

        cmd::varargxxx<util::pair_expression<int,int,':'>,double,double> eigenenergy(arg,"-eig","--eige","eigen energies of i:j part of the input hamiltonian within absolut e-range, i:j emin emax","",-1.0e+10,1.0e+10);

        cmd::varargxxx<util::pair_expression<int,int,':'>,double,double> scissor(arg,"-sc","--scissor","homo-lumo lev. shift of with dH and dL sel. range (cut at EF),i:j dH dL","",0.0,0.0);
        cmd::varargxxxx<util::pair_expression<int,int,':'>,double,double,double> scissormef(arg,"-scm","--scissorm","homo-lumo lev. shift with dH and dL of sel. range (cut at MEF),i:j dH dL MEF","",0.0,0.0,0.0);
        cmd::varargxxx<util::pair_expression<int,int,':'>,double,double> ascissor(arg,"-asc","--ascissor","homo-lumo lev. shift to H and L of sel. range (cut at EF),i:j H L","",0.0,0.0);
        cmd::varargxxxx<util::pair_expression<int,int,':'>,double,double,double> ascissormef(arg,"-ascm","--ascissorm","homo-lumo lev. shift to H and L of sel. range (cut at MEF),i:j H L MEF","",0.0,0.0,0.0);
//        cmd::varargxxx<util::pair_expression<int,int,':'>,double,double> scissorgap(arg,"-scg","--scissorg","homo-lumo level shift of sel. range (cut at fermi en),i:j Gap alpha,","",0.0,0.5);
//        cmd::varargxxxx<util::pair_expression<int,int,':'>,int,double,double> scissorn(arg,"-scn","--scissorn","homo-lumo lev. shift of sel.range (cut bw. lev. n-1 and n),i:j n H L","",0,0.0,0.0);
//        cmd::varargxxxx<util::pair_expression<int,int,':'>,double,double,double> scissoref(arg,"-scef","--scissoref","homo-lumo level shift of sel. range (cut at fermi en given here),i:j ef H L","",0.0,0.0,0.0);
//        cmd::varargxxxx<util::pair_expression<int,int,':'>,double,double,double> scissorgapef(arg,"-scgef","--scissorgef","homo-lumo level shift of sel. range (cut at fermi en given here),i:j Gap alpha ef","",0.0,0.5,0.0);
        
        cmd::switcharg<false> gollum(arg,"-g","--gollum-format", "hamiltonian file in gollum format");
        cmd::switcharg<false> phonon(arg,"-dyn","--dynmat-format", "dynamical matrix file as hamiltonian file");
        cmd::switcharg<false> squareformat(arg,"-sqmat","--sqmat-format", "square matrix format as hamiltonian file");

        cmd::switcharg<false> symmetrize(arg,"-s","--sym", "symetrize hamiltonian file before calculation");
        cmd::switcharg<false> adv(arg,"-a","--adv", "Advanced analysis and outputs");
//        cmd::varargxx<util::pair_expression<int,int,':'>,util::pair_expression<int,int,':'> > zeroing1(arg,"-z1","--zero1","zero hamiltonian elements between i:j and k:l","","");
//        cmd::varargxx<util::pair_expression<int,int,':'>,util::pair_expression<int,int,':'> > zeroing2(arg,"-z2","--zero2","zero hamiltonian elements between i:j and k:l","","");

        cmd::posarg<std::string> ham_fn(arg,"hamiltonain input file name", "");

        if(*help)
        {
                arg.print_help(std::cout,"transportab");
                std::cout << "examples:" << std::endl;
                std::cout << "transportab hamiltonian.dat -ha 3:4 -hb -4:-3 -n 100" << std::endl;
                std::cout << "transportab hamiltonian.dat -ha 3:4 -hb -4:-3 -n 100 -r 0.001" << std::endl;
                std::cout << "transportab -chain 0_4_0_4_0_0.1_-0.3_0.1_0.3_0.1_0_4_0_4_0 -ha 1:2 -hb -1:-2 -n 100 -scg 4:5 1.0 0.5"<< std::endl;
                std::cout << "version 0.99b" << std::endl;
                exit(EXIT_SUCCESS);
        }
        if(*info)
        {
                std::cout << "version 0.99b" << std::endl;
                std::cout << "Copyright (c) 2012 David Zsolt Manrique (david.zsolt.manrique@gmail.com)"  << std::endl;
                
                //return EXIT_SUCCESS;
                exit(EXIT_SUCCESS);
    
        }                                                                                                                
        
        
        //if(!(~ham_fn || ~chain)) throw std::runtime_error("Give hamiltonian file!");
        
        blockSH sh;

        if (~chain) interpret_chain(*chain,sh);
        else if(~ham_fn) 
        { 
    	    if(*gollum) load_gollum(*ham_fn,sh);
    	    else if (*phonon) load_dynamicalmatrix(*ham_fn,sh);
    	    else if (*squareformat) load_sqmat(*ham_fn,sh);
    	    else load(*ham_fn,sh);
	}
        //else if (~chain) interpret_chain(*chain,sh);
        else throw std::runtime_error("Give hamiltonian file!");
        
        if(*symmetrize) sh.symetrize();
        else if(sh.antisym_norm() > 1e-6) throw std::runtime_error("Hamiltonian or Overlap is not hermitian!"); /*check if gramm matrix e?*/

	if(~efermi) {sh.ef = *efermi;}
        
        //std::cout << sh.H << std::endl;
        //std::cout << sh.S << std::endl;



        blockSH sh_la,sh_lb;
        if(~la) { if(!*gollum) load(la.var2,sh_la); else load_gollum(la.var2,sh_la);  }
        if(~lb) { if(!*gollum) load(lb.var2,sh_lb); else load_gollum(lb.var2,sh_lb);  }

        //////////////////////////////////////////////////////////////////

        int ia,ja;
        if(~ha) { ia=std::get<0>(*ha); ja=std::get<1>(*ha); sh.reindex(ia,ja); }
        else if(~ah) { ia=std::get<0>(*ah); ja=std::get<1>(*ah); sh.reindex(ia,ja); }
//        else if(~fa) { ia=std::get<0>(fa.var1); ja=std::get<1>(fa.var1); sh.reindex(ia,ja); }
        else if(~la) { ia=std::get<0>(la.var1); ja=std::get<1>(la.var1); sh_la.reindex(ia,ja); }
        else if(~eigenenergy && !~nume);
        else throw std::runtime_error("Give -ha -ah range or -la range lead ham. file!");

        int ib,jb;
        if(~hb) { ib=std::get<0>(*hb); jb=std::get<1>(*hb); sh.reindex(ib,jb);  }
        else if(~bh) { ib=std::get<0>(*bh); jb=std::get<1>(*bh); sh.reindex(ib,jb);  }
//        else if(~fb) { ib=std::get<0>(fb.var1); jb=std::get<1>(fb.var1); sh.reindex(ib,jb); }
        else if(~lb) { ib=std::get<0>(lb.var1); jb=std::get<1>(lb.var1); sh_lb.reindex(ia,ja);  }
        else if(~eigenenergy && !~nume);
        else throw std::runtime_error("Give -hb -bh range or -lb range lead ham. file!");

        std::cout << "# --- transportab --- " << std::endl;
        std::cout << "#                                 " << std::endl;
        std::cout << "# HHH = .. - HA2 - HA1 - HS - HB1 - HB2 - .. "<< std::endl;
        std::cout << "#                                 " << std::endl;

        if(~eigenenergy)
        {
            int imm,jmm;
            imm=std::get<0>(eigenenergy.var1); jmm=std::get<1>(eigenenergy.var1); sh.reindex(imm,jmm);
	    
	    double emi = eigenenergy.var2;
	    double ema = eigenenergy.var3;
            std::cout << "# Eigen energies calculated of HHH[" << imm+1 << ":" << jmm+1 << "] within the (absolute)energy range [" << emi <<";" << ema << "] eV " << std::endl;
            
            SH shmm;
            sh.get_part(imm,jmm,shmm);

            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> evec;
            Eigen::VectorXcd eval;
            Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(shmm.S.inverse()*shmm.H);
            evec = eigen_solver.eigenvectors();
            eval = eigen_solver.eigenvalues();

            std::cout << "#eig:    ";
            for(int i = 0; i < eval.rows(); i++)
            if((eval(i).real()<=ema)&&(eval(i).real()>=emi))
            std::cout << std::right << std::setw(18+4) << std::setprecision(10) << eval(i).real() << ' '; 
            std::cout << std::endl; 
            std::cout << "#eig-ef: ";
            for(int i = 0; i < eval.rows(); i++)
            if((eval(i).real()<=ema)&&(eval(i).real()>=emi))
            std::cout << std::right << std::setw(18+4) << std::setprecision(10) << eval(i).real()-sh.ef << ' '; 
            std::cout << std::endl; 
            
            if (!(~nume)) std::exit(0);
        }


/* 
        if(~zeroing1)
        {
            int imm,jmm;
            imm=std::get<0>(zeroing1.var1); jmm=std::get<1>(zeroing1.var1); sh.reindex(imm,jmm);

            int kmm,lmm;
            kmm=std::get<0>(zeroing1.var2); lmm=std::get<1>(zeroing1.var2); sh.reindex(kmm,lmm);

            std::cout << "# HHH[" << imm+1 << ":" << jmm+1 << "," << kmm+1 << ":" << lmm+1 << "] = 0 (symmetrically)" << std::endl;
            
            SH shmm;
            sh.get_part(imm,kmm,jmm,lmm,shmm);
            shmm.S *= 0.0;
            shmm.H *= 0.0;
            sh.set_part(imm,kmm,jmm,lmm,shmm);

            sh.get_part(kmm,imm,lmm,jmm,shmm);
            shmm.S *= 0.0;
            shmm.H *= 0.0;
            sh.set_part(kmm,imm,lmm,jmm,shmm);

        }
        if(~zeroing2)
        {
            int imm,jmm;
            imm=std::get<0>(zeroing2.var1); jmm=std::get<1>(zeroing2.var1); sh.reindex(imm,jmm);

            int kmm,lmm;
            kmm=std::get<0>(zeroing2.var2); lmm=std::get<1>(zeroing2.var2); sh.reindex(kmm,lmm);

            std::cout << "# HHH[" << imm+1 << ":" << jmm+1 << "," << kmm+1 << ":" << lmm+1 << "] = 0 (symmetrically)" << std::endl;
            
            SH shmm;
            sh.get_part(imm,kmm,jmm,lmm,shmm);
            shmm.S *= 0.0;
            shmm.H *= 0.0;
            sh.set_part(imm,kmm,jmm,lmm,shmm);

            SH shmm2;
            sh.get_part(kmm,imm,lmm,jmm,shmm2);
            shmm2.S *= 0.0;
            shmm2.H *= 0.0;
            sh.set_part(kmm,imm,lmm,jmm,shmm2);

        }
*/
        if(~scissormef)
        {
            int imm,jmm;
            imm=std::get<0>(scissormef.var1); jmm=std::get<1>(scissormef.var1); sh.reindex(imm,jmm);

            std::cout << "# Scissor correction of HHH[" << imm+1 << ":" << jmm+1 << "] around MEF=" << scissormef.var4 <<" eV " << std::endl;
            
            SH shmm;
            sh.get_part(imm,jmm,shmm);

            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> evec;
            Eigen::VectorXcd eval;
            Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(shmm.S.inverse()*shmm.H);
            evec = eigen_solver.eigenvectors();
            eval = eigen_solver.eigenvalues();
            double homo = -1.0e+10;
            double lumo = 1.0e+10;
            for(int i = 0; i < eval.rows(); i++)
            if(eval(i).real()<=scissormef.var4) 
    	    { if(homo<eval(i).real()) homo=eval(i).real(); } 
    	    else { if(lumo>eval(i).real()) lumo=eval(i).real(); }

            std::cout << "# Levels of HHH[" << imm+1 << ":" << jmm+1 << "] :" << std::endl;
            std::cout << "#  MEF-EF = " << scissormef.var4-sh.ef <<" eV" << std::endl;
            std::cout << "#  HOMO-EF = " << homo-sh.ef <<" eV" << std::endl;
            std::cout << "#  LUMO-EF = " << lumo-sh.ef <<" eV" << std::endl;
            std::cout << "#  HLGAP = " << lumo-homo <<" eV" << std::endl;

            for(int i = 0; i < eval.rows(); i++)
            if( eval(i).real() <= scissormef.var4 ) eval(i)+=scissormef.var2; else eval(i)+=scissormef.var3; 

            std::cout << "# Corrected levels :" << std::endl;
            std::cout << "#  new HOMO-EF = " << homo+scissormef.var2-sh.ef <<" eV" << std::endl; 
            std::cout << "#  new LUMO-EF = " << lumo+scissormef.var3-sh.ef <<" eV" << std::endl;
            std::cout << "#  new HLGAP = " << lumo+scissormef.var3-(homo+scissormef.var2) <<" eV" << std::endl;
            std::cout << std::endl;

            shmm.H = shmm.S * evec * eval.asDiagonal() * evec.inverse();
            sh.set_part(imm,jmm,shmm);
        }

        if(~scissor)
        {
            int imm,jmm;
            imm=std::get<0>(scissor.var1); jmm=std::get<1>(scissor.var1); sh.reindex(imm,jmm);

            std::cout << "# Scissor correction of HHH[" << imm+1 << ":" << jmm+1 << "] around MEF=EF=" << sh.ef <<" eV " << std::endl;
            
            SH shmm;
            sh.get_part(imm,jmm,shmm);

            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> evec;
            Eigen::VectorXcd eval;
            Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(shmm.S.inverse()*shmm.H);
            evec = eigen_solver.eigenvectors();
            eval = eigen_solver.eigenvalues();
            double homo = -1.0e+10;
            double lumo = 1.0e+10;
            for(int i = 0; i < eval.rows(); i++)
            if(eval(i).real()<=sh.ef) 
    	    { if(homo<eval(i).real()) homo=eval(i).real(); } 
    	    else { if(lumo>eval(i).real()) lumo=eval(i).real(); }

            std::cout << "# Levels of HHH[" << imm+1 << ":" << jmm+1 << "] :" << std::endl;
            std::cout << "#  MEF-EF = " << sh.ef-sh.ef <<" eV" << std::endl;
            std::cout << "#  HOMO-EF = " << homo-sh.ef <<" eV" << std::endl;
            std::cout << "#  LUMO-EF = " << lumo-sh.ef <<" eV" << std::endl;
            std::cout << "#  HLGAP = " << lumo-homo <<" eV" << std::endl;

            for(int i = 0; i < eval.rows(); i++)
            if( eval(i).real() <= sh.ef ) eval(i)+=scissor.var2; else eval(i)+=scissor.var3; 

            std::cout << "# Corrected levels :" << std::endl;
            std::cout << "#  new HOMO-EF = " << homo+scissor.var2-sh.ef <<" eV" << std::endl; 
            std::cout << "#  new LUMO-EF = " << lumo+scissor.var3-sh.ef <<" eV" << std::endl;
            std::cout << "#  new HLGAP = " << lumo+scissor.var3-(homo+scissor.var2) <<" eV" << std::endl;
            std::cout << std::endl;

            shmm.H = shmm.S * evec * eval.asDiagonal() * evec.inverse();
            sh.set_part(imm,jmm,shmm);
        }
        if(~ascissormef)
        {
            int imm,jmm;
            imm=std::get<0>(ascissormef.var1); jmm=std::get<1>(ascissormef.var1); sh.reindex(imm,jmm);

            std::cout << "# Scissor correction of HHH[" << imm+1 << ":" << jmm+1 << "] around MEF=" << scissormef.var4 <<" eV " << std::endl;
            
            SH shmm;
            sh.get_part(imm,jmm,shmm);

            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> evec;
            Eigen::VectorXcd eval;
            Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(shmm.S.inverse()*shmm.H);
            evec = eigen_solver.eigenvectors();
            eval = eigen_solver.eigenvalues();
            double homo = -1.0e+10;
            double lumo = 1.0e+10;
            for(int i = 0; i < eval.rows(); i++)
            if(eval(i).real()<=ascissormef.var4) 
    	    { if(homo<eval(i).real()) homo=eval(i).real(); } 
    	    else { if(lumo>eval(i).real()) lumo=eval(i).real(); }

            std::cout << "# Levels of HHH[" << imm+1 << ":" << jmm+1 << "] :" << std::endl;
            std::cout << "#  MEF-EF = " << ascissormef.var4-sh.ef <<" eV" << std::endl;
            std::cout << "#  HOMO-EF = " << homo-sh.ef <<" eV" << std::endl;
            std::cout << "#  LUMO-EF = " << lumo-sh.ef <<" eV" << std::endl;
            std::cout << "#  HLGAP = " << lumo-homo <<" eV" << std::endl;

	    ascissormef.var2-=homo;
	    ascissormef.var3-=lumo;
	    
            for(int i = 0; i < eval.rows(); i++)
            if( eval(i).real() <= ascissormef.var4 ) eval(i)+=ascissormef.var2; else eval(i)+=ascissormef.var3; 

            std::cout << "# Corrected levels :" << std::endl;
            std::cout << "#  new HOMO-EF = " << homo+ascissormef.var2-sh.ef <<" eV" << std::endl; 
            std::cout << "#  new LUMO-EF = " << lumo+ascissormef.var3-sh.ef <<" eV" << std::endl;
            std::cout << "#  new HLGAP = " << lumo+ascissormef.var3-(homo+ascissormef.var2) <<" eV" << std::endl;
            std::cout << std::endl;

            shmm.H = shmm.S * evec * eval.asDiagonal() * evec.inverse();
            sh.set_part(imm,jmm,shmm);
        }
        if(~ascissor)
        {
            int imm,jmm;
            imm=std::get<0>(ascissor.var1); jmm=std::get<1>(ascissor.var1); sh.reindex(imm,jmm);

            std::cout << "# Scissor correction of HHH[" << imm+1 << ":" << jmm+1 << "] around MEF=EF=" << sh.ef <<" eV " << std::endl;
            
            SH shmm;
            sh.get_part(imm,jmm,shmm);

            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> evec;
            Eigen::VectorXcd eval;
            Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(shmm.S.inverse()*shmm.H);
            evec = eigen_solver.eigenvectors();
            eval = eigen_solver.eigenvalues();
            double homo = -1.0e+10;
            double lumo = 1.0e+10;
            for(int i = 0; i < eval.rows(); i++)
            if(eval(i).real()<=sh.ef) 
    	    { if(homo<eval(i).real()) homo=eval(i).real(); } 
    	    else { if(lumo>eval(i).real()) lumo=eval(i).real(); }

            std::cout << "# Levels of HHH[" << imm+1 << ":" << jmm+1 << "] :" << std::endl;
            std::cout << "#  MEF-EF = " << sh.ef-sh.ef <<" eV" << std::endl;
            std::cout << "#  HOMO-EF = " << homo-sh.ef <<" eV" << std::endl;
            std::cout << "#  LUMO-EF = " << lumo-sh.ef <<" eV" << std::endl;
            std::cout << "#  HLGAP = " << lumo-homo <<" eV" << std::endl;

	    ascissor.var2-=homo;
	    ascissor.var3-=lumo;

            for(int i = 0; i < eval.rows(); i++)
            if( eval(i).real() <= sh.ef ) eval(i)+=ascissor.var2; else eval(i)+=ascissor.var3; 

            std::cout << "# Corrected levels :" << std::endl;
            std::cout << "#  new HOMO-EF = " << homo+ascissor.var2-sh.ef <<" eV" << std::endl; 
            std::cout << "#  new LUMO-EF = " << lumo+ascissor.var3-sh.ef <<" eV" << std::endl;
            std::cout << "#  new HLGAP = " << lumo+ascissor.var3-(homo+ascissor.var2) <<" eV" << std::endl;
            std::cout << std::endl;

            shmm.H = shmm.S * evec * eval.asDiagonal() * evec.inverse();
            sh.set_part(imm,jmm,shmm);
        }
	
	//save("matrix.log",sh);        
	if(sh.antisym_norm() > 1e-6) throw std::runtime_error("Hamiltonian or Overlap is not hermitian after matrix tempering!"); /*check if gramm matrix e?*/        
	
        int im,jm;
        if(~la)
        {
            SH h11,h12,h21,h22; sh_la.get_part_2x2(ia,ja,h11,h12,h21,h22);

            std::cout << "# LA = .. - HLA1 - HLA2 - .. " << std::endl;
            std::cout << "# [ HLA1 HLA2 ] = LA[" << ia+1 << ":" << ja+1 << "]" << std::endl;
            std::cout << "# ||HLA1 - HLA2|| = " << (h11.H-h22.H).norm()/(h11.H.norm()+h22.H.norm()) << std::endl;  // make the norm looks ok
            std::cout << "# HA = HLA2 + de*SLA2 " << std::endl;



            int na = ja-ia+1; 
            if(na%2 != 0) throw std::runtime_error("Error: -la range has to be even (j-i+1 : even)");
            
            SH h00;
            sh.get_part(1,na/2,h00);

            double de = (h00.H - h22.H).diagonal().array().sum().real() / h22.S.diagonal().array().sum().real();

            std::cout << "# energy shift : " << de << std::endl;
            
            a.h = h22.H + de * h22.S;
            a.s = h22.S;
            a.v = h12.H + de * h12.S;
            a.t = h12.S;
            
            ia = 1;
            ja = na/2;
            im = ja+1;
        }
        else if(~ha)
        {
            SH h11,h12,h21,h22; sh.get_part_2x2(ia,ja,h11,h12,h21,h22);

            std::cout << "# ||HA1 - HA2|| = " << (h11.H-h22.H).norm()/(h11.H.norm()+h22.H.norm()) << std::endl;
            std::cout << "# [ HA2 HA1 ] = HHH[" << ia+1 << ":" << ja+1 << "]" << std::endl;
            std::cout << "# HA = HA2" << std::endl;

            a.h = h11.H;
            a.s = h11.S;
            a.v = h12.H;
            a.t = h12.S;
            
            ja = ja - (ja-ia+1)/2;
            im = ja + 1;
        }
        else if(~ah)
        {
            SH h11,h12,h21,h22; sh.get_part_2x2(ia,ja,h11,h12,h21,h22);

            std::cout << "# ||HA1 - HA2|| = " << (h11.H-h22.H).norm()/(h11.H.norm()+h22.H.norm()) << std::endl;
            std::cout << "# [ HA2 HA1 ] = HHH[" << ia+1 << ":" << ja+1 << "]" << std::endl;
            std::cout << "# HA = HA1" << std::endl;

            a.h = h22.H;
            a.s = h22.S;
            a.v = h12.H;
            a.t = h12.S;
            
            ia = ia + (ja-ia+1)/2;
            im = ja+1;
        }

        if(~lb)
        {
            SH h11,h12,h21,h22; sh_lb.get_part_2x2(ib,jb,h11,h12,h21,h22);

            std::cout << "# LB = .. - HLB1 - HLB2 - .. " << std::endl;
            std::cout << "# [ HLB1 HLB2 ] = LB[" << ib+1 << ":" << jb+1 << "]" << std::endl;
            std::cout << "# ||HLB1 - HLB2|| = " << (h11.H-h22.H).norm()/(h11.H.norm()+h22.H.norm()) << std::endl;
            std::cout << "# HB = HLB2 + de*SLB2" << std::endl;

            int nb = jb-ib+1; 
            if(nb%2 != 0) throw std::runtime_error("Error: -lb range has to be even (j-i+1 : even)");
            
            SH h00;
            sh.get_part(sh.n-nb/2,sh.n-1,h00);

            double de = (h00.H - h11.H).diagonal().array().sum().real() / h11.S.diagonal().array().sum().real();
            
            std::cout << "# energy shift : " << de << std::endl;

            b.h = h11.H + de * h11.S;
            b.s = h11.S;
            b.v = h12.H + de * h12.S;
            b.t = h12.S;

            ib = sh.n-nb/2;
            jb = sh.n-1;
            jm = ib-1;
        }
        else if(~hb)
        {
            SH h11,h12,h21,h22; sh.get_part_2x2(ib,jb,h11,h12,h21,h22);

            std::cout << "# ||HB1 - HB2|| = " << (h11.H-h22.H).norm()/(h11.H.norm()+h22.H.norm()) << std::endl;
            std::cout << "# [ HB1 HB2 ] = HHH[" << ib+1 << ":" << jb+1 << "]" << std::endl;
            std::cout << "# HB = HB2" << std::endl;

            b.h = h22.H;
            b.s = h22.S;
            b.v = h12.H;
            b.t = h12.S;
            
            ib = ib + (jb-ib+1)/2;
            jm = ib - 1;
        }
        else if(~bh)
        {
            SH h11,h12,h21,h22; 
            sh.get_part_2x2(ib,jb,h11,h12,h21,h22);

            std::cout << "# ||HB1 - HB2|| = " << (h11.H-h22.H).norm()/(h11.H.norm()+h22.H.norm()) << std::endl;
            std::cout << "# [ HB1 HB2 ] = HHH[" << ib+1 << ":" << jb+1 << "]" << std::endl;
            std::cout << "# HB = HB1" << std::endl;


            b.h = h11.H;
            b.s = h11.S;
            b.v = h12.H;
            b.t = h12.S;

            jb = jb - (jb-ib+1)/2;
            jm = ib-1;
        }
/*        else if(~fb)
        {
            SH h11,h12,h21,h22; sh.get_part(ib,jb,ib,jb,h11);
            
            std::cout << "# FAKE layer is used for HB1." << std::endl;
            std::cout << "# HB1 = HB2 = HHH[" << ib+1 << ":" << jb+1 << "]" << std::endl;
            std::cout << "# HB = HB1" << std::endl;
            
            h12 = h11;
            h12.H.setIdentity();
            h12.S.setZero();

            b.h = h11.H;
            b.s = h11.S;
            b.v = fb.var2*h12.H;
            b.t = h12.S;

            jm = ib-1;
        }
*/
        SH shm,shva,shvb;
        sh.get_part(im,jm,shm);
        sh.get_part(ia,ja,im,jm,shva);
        sh.get_part(im,jm,ib,jb,shvb);    

        std::cout << "# HS = HHH[" << im+1 << ":" << jm+1 << "]" << std::endl;
        
        m.h = shm.H;
        m.s = shm.S;
        
        m.va = shva.H;
        m.ta = shva.S;
        
        m.vb = shvb.H;
        m.tb = shvb.S;

        if(~smallimag)
        {

            std::complex<double> ii(0,1.0); 
            //Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> simag = m.h;
            //simag.setIdentity();
            double eps = *smallimag;

            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> evec;
            Eigen::VectorXcd eval;
            Eigen::ComplexEigenSolver<Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > > eigen_solver(m.s.inverse()*m.h);
            evec = eigen_solver.eigenvectors();
            eval = eigen_solver.eigenvalues();

            for(int i = 0; i < eval.rows(); i++) eval(i)+=ii * eps;

            m.h = m.s * evec * eval.asDiagonal() * evec.inverse();

            std::cout << "# eigenvalue(HS) = egeinvalue(HS) + i*"<<*smallimag<< std::endl;        
            
        }

        if(~randa)
        {
            std::cout << "# HA = HA + random * "<<*randa<< std::endl;        
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> r = a.v;
            r.setRandom();
            a.v += *randa * r.real().cast<std::complex<double> > ();
        }

        if(~randb)
        {
            std::cout << "# HB = HB + random * "<<*randb << std::endl;
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> r = b.v;
            r.setRandom();
            b.v += *randb * r.real().cast<std::complex<double> > ();
        }

        if(~randab)
        {
            std::cout << "# HA = HA + random * "<<*randab<< std::endl;        
            std::cout << "# HB = HB + random * "<<*randab << std::endl;
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> ra = a.v;
            ra.setRandom();
            a.v += *randab * ra.real().cast<std::complex<double> > ();
            Eigen::Matrix<std::complex<double>,Eigen::Dynamic,Eigen::Dynamic> rb = b.v;
            rb.setRandom();
            b.v += *randab * rb.real().cast<std::complex<double> > ();
        }

        if(*homogenize)
        {
            std::cout << "# HA = HB = (HA+HB)/2 " << std::endl;
            a.h += b.h; a.h *= 0.5; b.h = a.h;
            a.s += b.s; a.s *= 0.5; b.s = a.s;
            a.v += b.v; a.v *= 0.5; b.v = a.v;
            a.t += b.t; a.t *= 0.5; b.t = a.t;
        }

        std::cout << "# HA HS HB will be used!" << std::endl;

        std::cout << "#                                 " << std::endl;
        std::cout << "# --- end transportab --- " << std::endl;

        tol_p = *tol;
        emin_p = *emin;
        emax_p = *emax;
        ef_p = sh.ef;
        nume_p = *nume;
        sqrte_p = *sqrtesampling;
        advanced = *adv;
}


#ifdef DEBUG
void out_complex_matrix(Eigen::Matrix< std::complex<double> , Eigen::Dynamic, Eigen::Dynamic >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        //os.flags(std::ios::scientific);
        os.flags(std::ios::fixed);
        if(comment)
        {
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "Re Im ... Re Im" << ' '; 
                os << std::endl;
        }
        for(int i = 0; i < m.rows(); i++)
        {
	        for(int j = 0; j < m.cols(); j++)
                {
                        os << std::right << std::setw(w+4) << std::setprecision(p) << real(m(i,j)) << ' '; 
                        os << std::right << std::setw(w) << std::setprecision(p) << imag(m(i,j)) << ' '; 
                }
                os << std::endl;                        
        }                
        os << std::endl;                        
}
#endif

#ifdef DEBUG
void out_complex_vector(Eigen::Matrix< std::complex<double> , Eigen::Dynamic, 1 >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        os.flags(std::ios::scientific);
        if(comment)
        {
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "Re Im" << ' '; 
                os << std::endl;
        }
        for(int i = 0; i < m.rows(); i++)
        {
                os << std::right << std::setw(w+4) << std::setprecision(p) << real(m(i)) << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << imag(m(i)) << ' '; 
                os << std::endl;                        
        }                
        os << std::endl;                        
}
#endif

#ifdef DEBUG
void out_real_matrix(Eigen::Matrix< double , Eigen::Dynamic, Eigen::Dynamic >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        os.flags(std::ios::scientific);
        for(int i = 0; i < m.rows(); i++)
        {
	        for(int j = 0; j < m.cols(); j++)
                {
                        os << std::right << std::setw(w+4) << std::setprecision(p) << m(i,j) << ' '; 
                }
                os << std::endl;                        
        }                
        os << std::endl;                        
}
#endif


#ifdef DEBUG
void out_ssd(double e, Eigen::Matrix< std::complex<double> , Eigen::Dynamic, Eigen::Dynamic > & s,std::ostream & os,bool comment=true, int w = 12,int p = 5)
{
        os.flags(std::ios::scientific);
        if(comment)
        {
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "E-Ef [eV]" << ' '; 
                os << std::left << std::setw(w) << std::setprecision(p) << "maximum norm of S*adj(S)-I " << ' '; 
                os << std::endl;
                os << '#';
                os << std::left << std::setw(w-1) << std::setprecision(p) << "Re Im ... Re Im" << ' '; 
                os << std::endl;
        }
                
        os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << e << ' ';
        if( (s.rows() != 0) && (s.cols() != 0) )
        {
            Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic > m;
            m = s*s.adjoint();

            os << std::right << std::setw(w) << std::setprecision(p) << std::scientific << ((m-Eigen::Matrix< std::complex<double>, Eigen::Dynamic, Eigen::Dynamic >::Identity(m.rows(),m.cols())).array().abs()).maxCoeff() << ' ';
            os << std::endl;                
            for(int i = 0; i < m.rows(); i++)
            {
    	            for(int j = 0; j < m.cols(); j++)
                    {
                            os << std::right << std::setw(w+4) << std::setprecision(p) << std::fixed << real(m(i,j)) << ' '; 
                            os << std::right << std::setw(w) << std::setprecision(p) << std::fixed << imag(m(i,j)) << ' '; 
                    }
                    os << std::endl;                        
            }                
        }
        else os << std::endl;                        
        os << std::endl;                        
}
#endif

void out_complex_vector_as_line(const std::string &label,double energy,Eigen::Matrix< std::complex<double> , Eigen::Dynamic, 1 >  m,std::ostream & os,bool comment = true,int w = 18,int p = 10)
{
        os.flags(std::ios::scientific);
        os << '#' << label << ": ";
        os << std::right << std::setw(w) << std::setprecision(p) << energy << ' '; 
        for(int i = 0; i < m.rows(); i++)
        {
                os << std::right << std::setw(w+4) << std::setprecision(p) << real(m(i)) << ' '; 
                os << std::right << std::setw(w) << std::setprecision(p) << imag(m(i)) << ' '; 
                                        
        }                
        os << std::endl;                        
}

